The Two-Dimensional Phase-Unwrapping Problem
Two-dimensional phase-unwrapping is an important data-processing step in several important applications, including interferometric synthetic aperture radar (SAR), magnetic resonance imaging, optical interferometry, and adaptive optics. For an introduction, see D. C. Ghiglia and M. D. Pritt, “Two-Dimensional Phase Unwrapping: Theory, Procedures, and Software,” John Wiley & Sons, 1998.
The phase-unwrapping problem arises when a set of “phases,” which are continuous-valued variables, which can take on any real value, are only known modulo 2π radians, and it is desired to know the true absolute values of the phases. In the two-dimensional version of the phase-unwrapping problem, the phases are arranged in a two-dimensional, rectangular lattice. The rectangular lattice of all the phases is referred to herein as an “image” and one particular lattice point in the image will be referred to as a “pixel.”
In many applications, the phases serve as a proxy for some other physical quantity, such as time delay, and can be used to infer some desired information. For example, in the interferometric SAR application, the phases indicate a relative time difference for two radar signals to return from a target. Knowing the geometry of the radar system, one can use the phases to infer a terrain elevation map. However, because the phases are only measured modulo 2π radians, the phases have to be “unwrapped” to recover a sensible terrain elevation map.
Phases are often measured in units of radians, but it is also possible, and sometimes more convenient, to measure them in units of cycles. For example, if a particular phase has a value of 4.5 cycles, it could equivalently be said to have a value of 4.5(2π) radians. The convention of measuring phases in units of cycles is used herein. Therefore, an unwrapped phase has a value that can be any real number of cycles, while a wrapped phase always has a value in the interval [0.0, 1.0). Note that this standard mathematical notation means that 0.0 and any number greater than 0.0 but less than 1.0 is a possible value, but 1.0 is not.
FIG. 1 shows an example of three images 101-103 of 12 pixels arranged in a 3×4 rectangle. In FIG. 1, the circles represent pixels, and the numbers the pixel values. Of course, images of practical interest are much larger. The image 101 contains the absolute phases. In the image 102, all the phases are “wrapped” by adding or subtracting an integral number of cycles so that the phases are within the interval [0.0, 1.0).
Noise
In practical systems where the phase-unwrapping problem arises, one must also deal with the fact that the measurements of the phases are nearly always corrupted by noise. In the image 103, some of the wrapped phase measurements are corrupted by noise. This is indicated by darkened circles.
Criteria for a Good Phase-Unwrapping Method
A good phase-unwrapping system should both remove noise and unwrap the phases to recover an image that is as close as possible to the image of the true absolute phases. The goal is to take as input a noisy phase-wrapped image, such as image 201, and return as output an unwrapped image 202. A sensible criterion for selecting among the possible output solution images is to select that image which is most probable. Of course, to determine the probability of any solution image, one needs a model that assigns probabilities to the unwrapped phase image 202, given the noisy wrapped phase image 201.
In many applications, the images are extremely large. It is not very unusual for the noisy and wrapped phase image to contain more than 100,000,000 pixels. To be practical for such problems, a phase-unwrapping method must therefore be very fast, and the time to produce the absolute phase image should scale linearly or nearly linearly with the number of pixels.
To summarize then, a good phase-unwrapping method should use a well-justified probabilistic model of phase images, and it should efficiently find the most probable or nearly most probable absolute phase image according to that probabilistic model.
Factor graphs
Probabilistic models can be represented using a bi-partite graph that is called a “factor graph,” see F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor Graphs and the Sum-Product Procedure,” IEEE Transactions on Information Theory, vol. 47, pp. 498-519, February 2001. A factor graph contains two types of nodes, called “variable nodes” and “factor nodes.” Variable nodes are only connected to factor nodes and vice-versa.
The variable nodes represent the variables of interest in a system, while the factor nodes represent local functions of the variables. When a variable node is connected to a factor node, it means that the corresponding variable is an argument of the corresponding local function.
A probabilistic model is represented by a factor graph according to the equation
                              p          ⁡                      (                                          x                1                            ,                              x                2                            ,                              x                3                            ,              …              ⁢                                                          ,                              x                N                                      )                          =                              1            Z                    ⁢                                    ∏                              a                =                1                            M                        ⁢                                                  ⁢                                                            f                  a                                ⁡                                  (                                                            {                      x                      }                                        a                                    )                                            .                                                          (        1        )            
In equation (1), it is assumed that there are N variables that are denoted by x1, x2, . . . xN, and M local functions that are denoted by ƒ1, ƒ2, . . . , ƒM. The variables that are the arguments of the function ƒa are denoted by {x}a. A normalization constant Z ensures that the sum of the probabilities of all possible configurations is one.
An alternative but equivalent way to interpret a factor graph is in terms of “costs,” sometimes also called “energies,” rather than local functions. In this interpretation, the overall probability is given by
                                          p            ⁡                          (                                                x                  1                                ,                                  x                  2                                ,                                  x                  3                                ,                …                ⁢                                                                  ,                                  x                  N                                            )                                =                                    1              Z                        ⁢                          exp              ⁡                              (                                  -                                                            ∑                                              a                        =                        1                                            M                                        ⁢                                                                  c                        a                                            ⁡                                              (                                                                              {                            x                            }                                                    a                                                )                                                                                            )                                                    ,                            (        2        )            
where ca({x}a) is the cost for the variables {x}a to be in a given configuration, which is related to the local function ƒa({x}a) by an equationca({x}a)=−1n ƒa({x}a).
It is often more convenient and intuitive to describe factor graphs in terms of the “cost” formulation, because then all individual costs can be summed to obtain the total cost. A most probable configuration of the variables is the configuration with lowest total cost.
Prior Art Phase-Unwrapping Methods
A large number of prior-art methods for unwrapping noisy two-dimensional phase-wrapped images are known. However, few of the prior-art methods use a well-justified and explicit probabilistic model. Even fewer methods are capable of both unwrapping and de-noising phase images. The methods described herein are those that are relevant to the present invention.
The most similar prior-art method is the “ZπM” method described by J. M. B. Dias and J. M. N. Leitao, see J. M. B. Dias and J. M. N. Leitao, “The ZπM Procedure: A Method for Interferometric Image Reconstruction in SAR/SAS,” IEEE Transactions on Image Processing, vol. 11, pp. 408-422, April 2002. The ZπM method uses an iterative approach. In each iteration, the method uses a discrete optimization step followed by a continuous optimization step.
The probabilistic model used in the ZπM method uses a cost function that is quadratic in the difference between the absolute phases of two neighboring pixels. The cost function used in the ZπM method also depends sinusoidally on a difference between the input phase and the re-wrapped value of the obtained absolute phase.
Because of the form of the cost function used in the ZπM method, neither the discrete nor the continuous optimization steps are easy to perform. Therefore, Dias et al. describe complicated approximation techniques to perform the optimization steps.
A very important problem with the probabilistic model used in the ZπM method is that it does not explicitly constrain the phase derivatives to satisfy an “irrotational” constraint. This constraint is also known as the zero-curl constraint, see Dias et al. at p. 411. Dias et al. state that “this constraint is indirectly enforced by the prior [probability]: if the referred rotational is not zero at some point this implies the presence of large phase differences . . . which are penalized by the prior [probability].” However, this justification for dropping a very important constraint is not satisfying, because the phase difference cost is supposed to reflect a different aspect of the problem, i.e., the “smoothness” of the underlying phase image, rather than the inherent hard constraint represented by the zero-curl constraint.
Another prior-art method is called the “minimum cost flow” (MCF) method, see Constantini, M., “A novel phase unwrapping method based on network programming,” IEEE Transactions on Geoscience and Remote Sensing. 36, pp. 813-821, 1998, incorporated herein by reference. The MCF method suffers from the very serious drawback that it cannot correct any noise in the phase images.